Lattice Wess-Zumino model with Ginsparg- Wilson fermions: 
One-loop results and GPU benchmarks 

Chen Chenjj Eric Dzienkowskij^ and Joel Giedtlil 

Department of Physics, Applied Physics and Astronomy, 
Rensselaer Polytechnic Institute, 110 8th Street, Troy NY 12065 USA 

(Dated: July 30, 2010) 

Abstract 

We numerically evaluate the one-loop counterterms for the four-dimensional Wess-Zumino model 
p^ , formulated on the lattice using Ginsparg- Wilson fermions of the overlap (Neuberger) variety, to- 

gether with an auxiliary fermion (plus superpartners), such that a lattice version of U{1)r symmetry 
is exactly preserved in the limit of vanishing bare mass. We confirm previous findings by other 
P^ ' authors that at one loop there is no renormalization of the superpotential in the lattice theory, but 

that there is a mismatch in the wavefunction renormalization of the auxiliary field. We study the 
C^ \ range of the Dirac operator that results when the auxiliary fermion is integrated out, and show 

JL ! that localization does occur, but that it is less pronounced than the exponential localization of 

^ . the overlap operator. We also present preliminary simulation results for this model, and outline 

a strategy for nonperturbative improvement of the lattice supercurrent through measurements of 
CN ■ supersymmetry Ward identities. Related to this, some benchmarks for our graphics processing 

\^Q \ unit code are provided. Our simulation results find a nearly vanishing vacuum expectation value 

pg \ for the auxiliary field, consistent with approximate supersymmetry at weak coupling. 
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I. INTRODUCTION 

One would like to have a general method for studying strongly coupled supersymmetric 
field theories with lattice techniques. This is because nonperturbative dynamics play an 
important role in the theory of supersymmetry breaking and its transmission to the visible 
sector of particle physics. In this article we examine one such general method, which involves 
a fine-tuning of bare lattice parameters, after having restricted the number of counterterms 
using lattice symmetries. At the same time, we perform detailed numerical studies of a 
lattice formulation that was studied by several groups a few years ago [l|-l6|. We will highlight 
some interesting features of that model and present some new results regarding locality of 
the lattice Dirac operator and the degree of explicit supersymmetry breaking that occurs. 

Four- dimensional supersymmetric models on the lattice^ generically require fine-tuning 
of counterterms. This is to be contrasted with lower dimensional theories where lattice 
symmetries can eliminate the need for such fine-tuning; see [lO|] for further details. The 
one known four- dimensional exception is pure M = 1 super- Yang-Mills using Ginparg- 



Wilson fermions; the domain wall variety has been the subject of past 13|] and recent 
[qI, Il4|-|l8| simulations. Clearly we would like to go beyond pure M = 1 super- Yang-Mills. 
Recently it was proposed [l9| that an acceptable amount of fine-tuning could be efficiently 
performed using a mu 
Swendsen reweighting 



ticanonical Monte Carlo 20(| simulation together with Ferrenberg- 
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23| in a large class of theories; see also J9|. In this approach it is 
necessary to design the multicanonical reweighting function. We had suggested beginning 
at weak coupling on small lattices, using the one-loop perturbative counterterms as initial 
conditions to an iterative search approach. The present article reports numerical results 
for one-loop calculations in lattice perturbation theory that are designed to locate this 
starting point. Ironically we find that for the lattice theory that is studied here, the one- 
loop counterterms are entirely wavefunction renormalization, and that the logarithmically 
divergent parts match for the scalar and fermion. As a result, the initial condition for the 
iterative search is equivalent to starting with the tree-level action, since it is just a rescaling 
of the fields. This lends interest to our simulation results for the action with no fine-tuning, 
which we report here. 

The theory that we study is the four-dimensional Wess-Zumino model, formulated on the 



lattice with a variant of overlap (Neuberger) fermions 2J]. The goal of the formulation is 



;o impose the Majorana condition and simultaneously preserve the chiral U{1)r symmetry 
iMo] that is present in the continuum in the massless limit. As will be seen, preserving 
this symmetry greatly limits the number of counterterms that must be fine-tuned in order 
to obtain the supersymmetric continuum limit. In addition to overlap fermions, the lattice 
formulation has auxiliary fermions (plus superpartner fields) that couple to the overlap 
fermions through the Yukawa coupling, as in J25|. It is possible to integrate out the auxiliary 



^ For reviews with extensive references see |7Hl2l|. 



ferinions (and superpartner fields), and when one does this a nonanalytic Dirac operator 
results for the surviving fermionic field. Thus, as has been discussed originally in [l|, and 
at greater length in |2|, |5| the action is singular once auxiliary fields are integrated out. 
However, as we describe below, there is a sensible resolution of this singularity by taking 
the theory to "live" inside a finite box, with antiperiodic boundary conditions in the time 
direction for the fermions. The singularity of the Dirac operator that this resolves is related 
to nonpropagating modes in the infinite volume limit; the fact that these are nonpropagating 
shown in 6|. However, singularities in the Dirac operator raise the spectre of possible 



was 



nonlocalities in the continuum limit, as was found in gauge theories with the SLAG derivative 
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By analogy to that study, we have analyzed the continuum limit of the scalar self- 
energy and find that it is analytic in p^, so that there is no sign of nonlocality. We have also 



measured the degree of localization of the Dirac operator following the approach of 27| . We 
find that while there is localization, it is less pronounced than the exponential localization of 
the overlap operator. In the process of discussing our numerical perturbative results we are 
able to highlight the divergence structure of this theory, which turns out to be strictly wave 
function renormalization at one-loop. The wave function renormalization of the fermion and 
the physical scalar match at one loop in the continuum limit of the lattice expressions; but, 
the auxiliary scalar has a mismatched wave function counterterm. These findings appeared 
previously in [jj; thus we confirm those results. 

Having discussed the perturbative results, we then review correlation functions and renor- 
malization constants that must be measured in the simulations in order to fine-tune the 
theory. These involve the renormalized supercurrent (the current is renormalized because 
the symmetry is broken by the lattice regulator). We conclude with preliminary simulation 
results. In particular, we have developed all of the components for simulations on graphics 
processing units. Benchmarks that characterize the performance of our code are reported 
here. We measure one broken Ward-Takahashi identity in a simulation and find it to be 
very small at weak coupling. We argue that this is consistent with the nonrenormalization 
of the superpotential at one loop in lattice perturbation theory. 

II. DEFINITIONS 
A. Continuum 

The Euclidean continuum theory has action 

S=- f d'^x I ]-x^CMx + 0*n0 + F*F + F*(m*0* + g^cj)*^) + F(m0 + ^02)1, 

M = ^+{m + 2g<p)P+ + (m* + 2g*(p*)P_. (2.1) 



Our conventions will be {i = 1, 2, 3): 

P± = ^{l±l5), C = 7o72. (2.2) 

It can be checked that the action is invariant under the supersymmetry transformations 

6,(P = y2e^CP+x, Se(p* = y2e^CP_x, 

6,x = -v^P+(# + F)e- V2P_(#* + F*)e, 

5,F = V2e^C^P+x, 5eF* = V2e^C^P_x (2-3) 

Note also that in momentum space 

-0*n0^|0(p)lV (2.4) 

so that for a well defined partition function we must take the negative sign in the exponent: 

Z = f[d(P dcj)* dF dF* dx] e-^ (2.5) 

On the other hand, integration over the auxiliary field F is not well-defined (as is usual in 
Euclidean formulations of supersymmetric theories). In particular, we have after integrating 
out the fermions and completing the square, 

Z = f[d(f) rf</)*]e-^(*)pfCM(0) f[dF dF*] exp f d^x \F* + m(j) + g(f)'^\'^, 

5(0) = f d'^x (-0*n0+|m0 + ^02|2) (2.6) 

where Pf CM is the Pfafiian of the matrix CM. Under shift of integration variables 

F ^ F' = F + m*(j)* + g*(j)*^, F* ^ F'* = F* + mcj) + gcj)'^ (2.7) 

Z does not converge, due to the auxiliary field factor j{dF' dF'*] exp J d'^x |P'p. This fact 
was noted, for example, in [l| where they simply divide by an infinite constant — the partition 
function of the (7 = theory — to cancel the infinity. In a lattice simulation we do not have 
this luxury, and we must decide on the proper way to deal with the integration over F, F*. 
Of course in the Minkowski space formulation, due to the presence of an additional factor 
—2, the integral over the auxiliary fields F', F'* can be computed by analytic continuation. 
One discards the overall constant into the normalization of Z, and the net result is that one 
imposes the equations of motion: 

^^ ^^ 0. (2.8) 



5F{x) 5F*{x) 



From this we conclude that in Euchdean space the integral fl2.6p should be understood in 
a formal sense; it is an instruction to impose fl2.8p . This is equivalent to completing the 
square, as in fl2.6p . shifting the integration variable and then ignoring the (infinite) constant 
factor that is generated. Thus all simulations are performed with the action in the form with 
the auxiliary fields removed. In this form the supersymmetry transformations are nonlinear 
and the supersymmetry algebra only closes on-shell. 



B. Lattice 



We next discuss the lattice action, which is a special case of the formulations of [l|, l2| ; 



write the lattice action in forms given in 
must be introduced. 



we 



6|. For this, a few lattice derivative operators 



1 1 

A = I- aDw, Dw = ipt,{dl + <9^) + -j^ad^d^ 



D2 
D 



1 + l^X^ ] (^^^) 



-1/2 



D^ + D2 = -{l-A{A^Ay'/^) 



(2.9) 



where 9^ and d* are the forward and backward difference operators respectively. Then the 
lattice action is Isl: 



S = -a^J2{ \x^'^^X + <P*Dl(p + F*F + FD^cf) + F^D^cf)* 

X ^ 

--X^CX - - (J^$ + J^*<l>*) 
a a 



F*{m*(j)* + g*(j)*^) + F{m(f) + 



Here, the tilded fields are the linear combinations 



(2.10) 



+ $, x = X + X, F = F + J^ 



(2.11) 



The fields $, X, T and their conjugates are auxiliary fields introduced to allow for a lattice 
realization of the chiral ?7(1)_r symmetry in the m — )■ limit: 



(50 = 


— 3ia0 + ia 


[{.-lo.y 


-f-] 


H 


(5$ = 


-3i«$ + i^a (^20 + i^*) 




5F = 


3iaF + ia 


f a \ 
(l--fl.)F- 





+ zq; J-" 



(5J^ = 3iaJ^ + i-a {D2F + Dl(p*) 



(2.12) 
(2.13) 



(2.14) 



which takes a particularly simple form on the tilded variables: 

6x = ict'-y^x^ ^4> = —'^ic(4>, SF = AiaF 
The supersymmetry transformations of the untilded fields will be taken as 

(5,0 = y2e^CP+x, '5e0* = V2t^CP^x, 
5ap = -v^(P+(Di0 + F)e)p - v^(P„(Di0* + F*)e)p, 
6,F = V2e^CDiP+x, 5,F* = V2e^CDiP_x 
and for the auxiliary multiplet 

5,^ = V2e^CP+X, 6,^* = V2e^CP-X, 

5,Xp = -v^(P+(Di$ + ^)e)^ - v^(P_(Di$* + r)e)p, 

5^jr = y/2e^CDiP+X, 6,J^* = V2e^CDiP_X 

Of course this is not a symmetry of the lattice action for g ^ 0. Here we are just choosing the 
form of the tranformation that will be used to generate broken Ward-Takahashi identities 
on the lattice. The supersymmetry transformations of the tilded fields are 

6,4> = v^e^CP+x, 5e0* = V2e^CP^x, 

6,XP = -V2{P4Di4> + F)e)[s - V2(P_(Di<^* + P*)e)^, 

6,F = y/2e^CDiP+x, S,F* = V2e^CDiP^x (2.16) 



(2.15) 



As noted in [5|, we can integrate out the auxiliary fields X, <l>, J^, treating the tilded fields 
as constant, to obtain the lattice action:^ 

S = -a'Y,{lfCMx-l4>*D,4> + F*{l-lD2)-'F 

X ^ 

+F*{m*4>* + g*4>*^) + F{m4) + g^^)\. (2.17) 



^ Integrating out an auxiliary fermion to obtain the fermionic part of this action was previously noted 
in ^]. There it was noted that this relates the fermionic action to the one of [25| by a singular field 
transformation. This singularity will be discussed more in Section [IIII below. 



This is the lattice action Eq. (2.14) of [^ with a notation that interchanges Di ^ -D2, which 
is equivalent to Eq. (2.22) of |2| for the k = case, using the identities^ 

T, = j,{l-^D), Tl = l-^D,, D^D = ^D,. (2.18) 

The fermion matrix is: 

M = ]p + mP+ + m*P_ + 2g4>P+ + 2g*4)*P^, ]p = {1 - ^Diy^Di (2.19) 

This way of writing the Dirac operator can be related to the one that appears in [2,] by the 
identity: 

(1 - ^D,)-'D, = (1 - ^D)-'D (2.20) 

Furthermore we can integrate out the auxiliary fields F, F* to obtain the action 

^ = a'J2{-lx^CMx + l4>*D24>+{m*4>* + 9*4>*'){l-^D2){m4> + g4>')y2.2^^ 

This is the action that is used in our numerical simulations. 

When fine-tuning of the lattice action is performed, we must invoke the most general 
lattice action consistent with symmetries. Since we perform our simulations at m ^ 0, this 
is just the action with all dimension < 4 operators built out of the physical fields, and x- 
We write it here for reference: 

•^ = «' E i - k'^C'(^ + miP+ + mlP,)x + -^*D24> 
^ — ^ 2 a 

X ^ 

-X^C(yi0P+ + yl^*P-)x - fC{y,~^P^ + yfrP+)x\ (2.22) 



A term linear in cf) has been eliminated by the redefinition </> — )■ + c with c a constant. 
The parameters m\ and Ai are real and all other parameters are complex. Whereas in 
the supersymmetric theory there are four real parameters, in the most general theory we 
have eighteen real parameters to adjust. Holding mi and yi fixed, we have fourteen real 
parameters that must be adjusted to obtain the supersymmetric limit. The counting can 
be alleviated somewhat by imposing CP invariance, so that all parameters can be assumed 
real. Then we have a total of ten parameters. Holding two fixed, we must tune the other 
eight to achieve the supersymmetric limit. Conducting a fine-tuning in an eight-dimensional 
parameter space is a daunting task. 



We thank A. Feo for explaining this point and providing us with a derivation of these relations. 



On the other hand in the hmit mi — )■ we can impose the t^(l)i? symmetry f l2.13p . This 
restricts the action to 

-fC{yi^P+ + yl^*P^)x] (2.23) 



If we hold yi fixed, then only m^ and Ai must be fine-tuned. Conducting a search in a two- 
dimensional parameter space, with both coming from bosonic terms, is manageable. The 
difficult part is that we must extrapolate to the massless fermion limit. Another potential 
problem is that we impose antiperiodic boundary conditions for the fermion in the time 
direction, but must impose periodic boundary conditions for the scalar in order for the 
action to be single-valued on the circle in the time direction. This breaks supersymmetry 
explicitly by boundary conditions. At finite mass this should be an effect that can be made 
arbitrarily small by taking the large volume limit. However at vanishing mass, there will be 
long distance modes that will "feel" the breaking due to boundary conditions.^ Thus it is 
important that we take T ^ 1/ma as m is sent to zero, where T is the number of sites in 
the time direction. 

As we will see, the one- loop behavior of the theory ( I2.17P closely follows that of the 
continuum, so that no new operators are generated at this order. Thus at this level of 
approximation, a fine-tuning of the general lattice action f l2.22p is not needed. For reasons 
that will be clearer once we have presented our perturbative results, it is of interest to 
study the original lattice action fl2.2ip in our simulations, without any fine-tuning. By 
measuring the degree of supersymmetry breaking through supersymmetry Ward identities 
that are conserved in the continuum, we gain information about the higher orders and 
nonperturbative aspects of the lattice theory. 

III. MODES, FEYNMAN RULES AND TILDE/UNTILDE EQUIVALENCE 

A. Mode analysis 

As to the the nonlocal operator (1 — ^aD2)^^ that appears in the action ( I2.17p . note that 
for smooth field configurations aD2 ~ a^D. Nevertheless, one may rightly worry about the 
effect of modes for which D2 ~ l/a. The spectrum of 1 — ^aD2 can easily be calculated in 
momentum space: 

1 1 1 l-2E„sin^(pW2) 

2 2 2 /[1-2V sin2(p,a/2)]2 + V sin2(p,a) 



We thank G. Bergner for raising this point. 



One sees that there are zeros at p^a equal to all of the would-be doublers: 

(ttAM), ( 7r,7r,0,0 ), ( tt, tt, tt, ), (tt, tt, tt, tt), (3.2) 

where the underline indicates all possible permutations. Thus the auxiliary field kinetic 
term appearing in f l2.17p is singular, as was discussed in |l|, and more at length in (2|, |5|, |6[. 
For the fermions we must also consider 



-^« ^Eu7Msin(P;. 



a] 



DM = "^ ^ — (3.3) 

[1 - 2 E, sin2(p^a/2)]2 + ^^ sin\p^a) 

Thus Diij)) = for the would-be doublers (13.21) and the fermion operator (1 — ^aD2)^^Di 
is indeterminate. In |6|] it was shown that in a limiting process of approaching these points 
in momentum space, (1 — ^aD2)^^Di is divergent. Thus the would-be doublers are actually 
nonpropagating. Note that 

(1 - ^aD2)-'Di^5 + 75(1 - ^aD2)-'Di = (3.4) 

which is how the lattice formulation manages to preserve U{1)r symmetry. This is not 



in conflict with the Nielsen- Ninomiya theorem |28l . |29| | precisely because the operator (1 — 
^aD2)~^Di is unbounded and therefore nonanalytic. One should then worry about locality 
since a nonanalytic Dirac operator might violate this property, as was true of the SLAG 
Dirac operator [26|. We will investigate this below. We will find numerically that it does 
show localization, though not as pronounced as the exponential localization of the overlap 
operator. In fact it has a long tail that is a cause of some concern. However our perturbative 
results in a later section show no sign of nonanalytic sickness in the scalar self-energy as a 



function of p^, which would be the analogue of the results of 26 1. 

For the fermion, in our simulations, we address the difficulty of the indeterminate operator 
by taking the fermion field to be antiperiodic in the time direction and restrict considerations 
to finite lattices L^ x T with T even. Taking T even is necessary in order to avoid other 
zeros of (1 — ^aD2). This resolves the indeterminacies for all finite L,T with T even. One 
then defines the infinite spacetime volume theory rigorously as the limiting value L, T — )■ 00. 
Unfortunately, the maximum eigenvalue of the fermion operator (1 — ^aD2)^^Di diverges 
as T — )■ 00, so that the fermion matrix becomes poorly conditioned at very large system 
size. For the T = 32 and T = 64 lattices that we study this does not prove to be a serious 
problem. 

For the auxiliary field F, in our simulations, we formally integrate it out, leading to the 
action (I2.2ip . For the modes with (1 — ^D2){p) = the scalar potential vanishes. However, 
the kinetic term -0*1^20 is nonvanishing on this modes, with -1^2 (p) = 4/a^. These modes 
are therefore suppressed in the path integral and are part of the cutoff theory. 

In perturbation theory, one can arrange for the modes at p^a of the would-be doublers 
(13. 2p to be nonpropagating, as ought to be based on the results of [6| . This nonpropagating 



feature is incorporated into the Feynman diagram rules of [5], as we now explain. First we 
note that 

-1 



-Di + ( 1 - -D2 ) {m*P+ + mP. 



(l - ^Z^s) Di + {mP+ + mTP.) 



D, + \m\'(l-^D,) 



using the Ginsparg- Wilson relation, which in terms of Di and D2 takes the form 

2 



Di - D'^ = — L>2 

a 



Then we find that the fermion propagator can be written 

-Di + (1 - f Z^a) {m*P+ + mP.) 



{x{x)r {y))c = a 



lD2+\m\^{l-lD,) 



{x,y) 



(3.5) 



(3.6) 



(3.7) 



This propagator vanishes for the would-be doublers (13.21) . It avoids the alternative form that 
comes from a straightforward inversion of the free Dirac operator Mq = (l — I-D2) -Di + 
mP+ + m*P_: 



{xix)x' iy))c = -a- 



(1-1^2) ^Di-{m*P+ + mP_] 



{x,y) 



(3.8) 



which is indeterminate for the would-be doublers. In our perturbative analysis of Section 
[IV] below, we use the propagator (13. 7p . This allows us to use periodic boundary conditions 
for the one-loop calculations. 



B. Locality 



Following the numerical approach of 271] we can investigate the locality of the operator 
p = [1 — f-D2) -Di. In what follows, we impose the antiperiodic in time boundary con- 
ditions on the fermions. We introduce a unit point source rj at the origin a; = and then 
compute 



ip^x) 



l-^D,y'D,rj 



[X 



The "taxi-driver distance'' 



\x\ 



|i = ^i\XfM\ or \L^ - Xf,\) 



(3.9) 



(3.10) 



to the origin is defined; the shortest length is selected in each "or" statement, where L^/a 
is the number of lattice sites in the /x direction. One then computes the norm [[-^(x)!! in 
spinor index space at site x and obtains the function 



/(r) = max{||-?/'(a;) 



|a;||i = r} 



(3.11) 
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FIG. 1: Probe of the range of the Jf) operator, Eq. (j3.1ip . The value of / at r = is zero and is 
not shown in this plot. Even r are consistently lower than odd r. 



We show this function for a 16^ lattice in Fig. [H What one sees is that there is a long tail on 
the operator. In contrast to the overlap operator, we do not find exponential localization; the 
localization that does occur is less pronounced and we view it as an open question whether 
or not this is harmful. 



C. Tilde/untilde equivalence 

Suppose we wish to compute correlation functions of the untilded fields. For this purpose 
we introduce sources that couple only to them: 



S. 



^\,^CP^,.r^CP.,.,;..,,-,U;.,,-] 



(3.12) 



Next we integrate out the auxiliary fields X, $, T . The corresponding equations of motion, 
determined by variation of X, $, J-" holding x, 0, F constant, involve the sources. These can 
then be solved for the auxiliary fields X, $, J-": 



X = -^{l-^Dy\Dx + P+V + P-r] 



2 V 2 



$ = 

2 V 2 



a / a \ -1 

k* (1--L)2) J 

2 4 V 2 / 



(3.13) 
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Thus when these are used to ehminate X, $, J-", quadratic terms in the sources are generated. 
We obtain 

+ ^X^ (mP+ + m*P_ + 2(?0P+ + 2(7*0*P_) x 

+F*{m*(j)* + c/*0*^) + F(m0 + ^0^) 

+^(r]^n + r^P-)C (l - ^/^)"' {P^V + i^-r) 

+ (77^P+ + r^P_)C (l - ^Dy\ + k{F* - ^20+ ^ J) 

+r (P - P20* + ^ J*) + J^* + J*0 - ^A;*P2A; 

+^F[i-^D,y\j + ^F*[i-^D,y\r-^ijj + .r.r) 

+^J*(^l-2aD, + ^{Dl + Dl)^ (l-^D.y'j^ (3.14) 

As far as the source terms for the elementary fields x ^i-nd are concerned, working with 
the tilded fields is the same as working with the untilded fields up to 0{a) corrections. Note 
however that the source terms involving k, k* have unusual terms that do not vanish in the 
continuum limit. The implication is that correlation functions for P, P* differ from those of 
P, P* in a way that is guaranteed not to vanish in the continuum limit. 

To see what are the implications for the physical fields, we set k = k* = and integrate 
out the auxiliary fields P, P*: 

S + Ssrc = -«' E { ^^^^ (l - i^) ' ^^ - l^*^^^ 
X ^ 

+ ix^ (mP+ + m*P_ + 2g^P+ + 2(?*0*P_) x 
+ K0* + g*^*') (l - ^P^) (m0 + g4>') 
+ ^{r^^P^ + r^P_)C (l - ^P)"' (P+r^ + P^r) 
+ [rfp^ + r^P_)c{l-'^Dy\ 

+ [^* (0 - ^(^0 + ^0')) + h.c] - j( J' + r') 



0- ,-. I ^ a . ^2 , n2 



2J*(^P2--Pf + /^2')J (1--/^2J J| (3.15) 

To interpret this result, suppose instead we had introduced sources J, J*,fi,f for the tilded 
fields 0*,0,P+X,^-X: 

Ssrc = -a' Yl {f^P+X + f^CP-X + J*4> + ^0*1 (3.16) 

12 



From fl3.15p what we see is that 

oJ{x) SJ{x) 

and hkewise for the other sources J*, fj, f. This is another way of saying that the correlation 
functions for the tilded fields are equal to the correlation functions of the untilded fields, up 
to 0{a) corrections. Of course if there were 0{l/a) or worse divergences in the correlation 
functions, then the difference between the two sets of correlation functions does not vanish in 
the continuum limit. However, what we will find below is that at one loop the divergences are 
only In a, hence the tilded and untilded correlation functions become equal in the continuum 
limit. If 0{l/a) or worse divergences appear at higher orders, then the two formulations are 
not equivalent in the continuum limit. In that case one is free to choose one or the other as 
the subject for fine-tuning in order to approach the continuum theory. We will choose the 
tilded formulation in our work. 

D. Naive continuum limit 

For reference, we state the basic properties of the lattice theory the guarantee that the 
correct continuum limit is achieved classically. Two useful identities involving the lattice 
derivative operators are 

2 
lim Di = ^, lim -D2 = -D (3.18) 

a— >-0 a-^0 (X 

Then in the naive continuum limit 

(i-|^2)-^^(i + jn)-^^i (3.19) 

Thus we see that for the kinetic terms 



a^ 



J2 {lx^C{l - ^D,)-'D,x - 14>*D,4> + F*(l - ^aD^r'F 



^ ^d^x |^x^C^X + 0*n0 + ^*^} (3.20) 

It is interesting that the continuum limit of the action written in terms of the tilded fields 
is correct. This is a further indication that we can just as well treat them as the "physical" 
variables and work entirely in terms of the action 02.17P . 



IV. ONE-LOOP CALCULATION 



In pj an identical calculation is performed, in that they also compute one-loop cor- 
rections to the propagators and proper vertices. The action that they use, described by 
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their Eq. (2.17), is different from tlie one employed here, but is related by a singular field 
transformation [2|. We have confirmed the results of [l|, working in this alternative but per- 
turbatively equivalent formulation. Our results also lead to conclusions identical to those 
of |5|, who take a different approach to studying the same action as the one we use, our 
Eq. f l2.17p : they study the counterterms that must be added to the action in order for the 
restoration of the supersymmetry Ward-Takahashi identity to occur. 

A. Definitions 

Throughout our presentation, we make use of 

P^{k) = a-i sin(A;^a), Q^{k) = 2a-^ sm{k^a/2), (4.1) 

periodic functions that reduce to momentum in the naive continuum limit. It is also useful 
to define 

M(A;) = l-2^sin2(A;^a/2) = l-— g2(A;), s{k) = ^/P^{k)a^ + M^{k) (4.2) 
which has the property lim(j_>o s{}z) = 1. 



B. Analytic results 

The only nonvanishing one-loop scalar two-point function is the one that gives a one-loop 
correction to the nonholomorphic term in the effective action, 

d'x d'y 4>''ix)Gj!^{x,y)4>{y) (4.3) 

where G{x,y) is the scalar propagator. It is obtained from the sum of the two diagrams 
shown in Fig. [2], yielding: 

2, ,2 r^^ d'k N{k,p) ,,,, , n{k,p) 



P'^{k)a? 
n{k,p) = a^P{p + k)-P{k), d{k,p) = s{k)s{k + p), t{k) = ;./. (4.5) 

s^(/e) 

Dik^p) = rik)rik + p), Kfc) = 2(l-^)+^(l + ^) (4-6) 

The result vanishes at p = 0, because N{k, 0) = 0, so no mass counterterm arises from this 
diagram. 
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It is interesting that a term in the effective action of the form 

d\ dS mGlUx,y)^{y) (4.7) 



is not generated at one-loop, since it is allowed by the symmetries of the lattice action if 
m 7^ and (? 7^ 0. At one-loop there is an exact cancellation between the scalar and fermion 
loops for all external momentum p. 

The diagram that corrects the fermion propagator, Fig. [31 is given by: 

i-^M^ I ?^P ^(^ + ^) (4 8) 

^ ^'^' J (27r)4 - r{k)r{p + k)s{p + k) ^ ' 

The diagram Fig. H] gives the correction to the F*F term and is given by: 

^"^^'^' y {2T,Yr{k)r{p + k) ^^-^^ 

The 3-point diagram shown in Fig. [5] evaluates to: 

-im\g\''g* f a^ d^k fjp' + k) {1 + M{p + k)/s{p + k)) ^^ 

2a J (27r)4 r{p' + k)r{p + k)r{k)s{p' + k) " ^ ' ' 

Note that when the external momenta are set to p = p' = 0, this expression vanishes, due 
to an odd integrand, f{—k) = —f{k). Thus the Yukawa coupling receives no corrections in 
the zero-momentum subtraction scheme. 

The Z — 1 counterterms (cf. Eq. (14.161) below) are computed in the usual way, from the 
sum of amputated one-particle irreducible self energy diagrams S(p). For instance in the 
case of the scalar 

pZ _|_ ^ri2^ pZ _|_ jj^Z pZ _j_ jy^Z pZ _|_ .jj^Z pZ _j_ jy^Z 

where ■ ■ ■ represents higher order terms. Since there is no mass renormalization at one loop 
we have S(0) = and S(p) = 1(7^22 p^ + O^p"^). So, to have cancellation of E2 at p = we 
require Z^ — 1 = |5'pS2- 

C. Numerical results 

We have computed the loop integrals in two ways: working with discrete sums at finite 
system size L, and with numerical integration for infinite system size L — )■ 00. We set 
the scale using the bare mass m. In the discrete sums, we use mL = 8; the number of 
lattice sites in any direction is A^ = L/a = mL/ma. Equivalently, ma = rriL/N determines 
the lattice spacing in units of m. We have also computed with mL = 16 and find that 
the difference is only in the third or fourth significant figure. For the scalar we compute 
Z^ — 1 = limp_^o S(/i(p)/p^) with S^(p) the one-loop self-energy Fig. [2l since S(p) vanishes 
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FIG. 2: Scalar 2-point function that contributes to self-energy. 




FIG. 3: Fermion 2-point function that contributes to self-energy. 

at p = 0, as noted above. The results are given in Table HI For the fermion we compute 
Z^ — 1 = limp_j.o S^(|^)/|^ since ^x(^) vanishes at ^ = 0. The results are given in Table HTl 
For the auxiliary field, we compute Zp — 1 = Si?(0). The results are given in Table HTTl 
We have fit the L = oo numerical data with ma < 0.125 to 

f{ma) = Co In(ma) + Cima + C2(ma)^ (4-12) 

Giving the data points equal weight, the fit for the scalar self energy is 

Co = 0.00604(7), ci = 0.024(20), Ca = -0.15(17) (4.13) 

while for the fermion the fit is 

Co = 0.00597(5), ci = 0.061(15), Ca = -0.39(12) (4.14) 
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FIG. 4: Auxiliary 2-point function that contributes to self-energy. 




FIG. 5: 3-point function. 

Thus we see that at L — )■ oo the log divergences of the scalar Z^ — 1 and the fermion Z-^ — 1 
match. For the auxiliary field we obtain instead 



Co = -0.0261(5), ci = 0.87(11), Ca = -3.9(1.0) 



(4.15) 



so that its Zp — 1 does not match the scalar and fermion. This is consistent with the results 
found in |5|. 

D. Renormalization 



We can absorb all the logarithmic divergence by rescaling the fields, so that the action is 
modified from (I2.17P and (I2.19p to read 



S = —a 



^ J2 {l^^X^CMx - Z^4>*D24> + ZpF*{l - ^D,)-'F 



+ JZFF*(m* 



^ I *^ ^/ ^J^^ ^ ^ 



ZFF{m^JZ^(t) + gZ^(t) ) >, 



M = p + {m + 2g^Z^<P)P+ + (m* + 2gyZ^<P*)P^ 
^ = (1 - ]^aD2)-^Di 
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(4.16) 



ma 


{Z-l){oo)/\g\^ 


{Z-l){8)/\g\^ 


0.5 


-0.00589 


-0.006355 


0.25 


-0.00884 


-0.009532 


0.125 


-0.01205 


-0.013358 


0.0625 


-0.01573 


-0.017541 


0.03125 


-0.01986 


-0.021837 


0.015625 


-0.02573 


-0.026205 


0.0078125 


-0.02956 


-0.030583 


0.00390625 


-0.03370 


— 


0.001953125 


-0.03671 


— 



TABLE I: Self-energy for the scalar, from Fig. [2] evaluated at p = (0, 0, 0, 2tt/L) for L = oo and 
mL = 8. 



ma 


{Z-l){oo)/\g\^ 


{Z-l){8)/\g\^ 


0.5 


-0.00385 


-0.004027 


0.25 


-0.00715 


-0.006973 


0.125 


-0.01085 


-0.011406 


0.0625 


-0.01463 


-0.015632 


0.03125 


-0.01893 


-0.019965 


0.015625 


-0.02319 


-0.024336 


0.007813 


-0.02910 


-0.028720 


0.00390625 


-0.03311 


— 


0.001953125 


-0.03697 


— 



TABLE II: One-loop counterterm Z — \ for the fermion, evaluated from slope in self-energy in 
p — >• limit, for L = oo and m,L = 8. Note that the lattice spacing is measured in units of 1/m. 

It is interesting that when one eliminates F, F* by their equations of motion, the constant 
Zp disappears entirely: 



S 



Y.{\ 



— (m 



Z^x'CMx-Z^^*D2^ 



+ g*Z^<P*'){l--D2){m 



gZ^cj)^ 



(4.17) 



The discrepancy between Zp and Z^^ = Z^ is irrelevant to the physical theory, which only 
contains and x- Iii terms of the on-shell formulation, the one-loop renormalization of the 
lattice theory, in the a — )■ limit, exactly mirrors what occurs in the continuum theory. 
The mismatch Zp ^ Z^ = Z^ will have effects at two loops, where one-loop corrected 
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ma 


{Z-l){oo)/\g\^ 


(Z-l)(8)/|5p 


0.5 


0.100785 


0.082542 


0.25 


0.0994727 


0.088658 


0.125 


0.103388 


0.095266 


0.0833333 


0.106955 


0.09916 


0.0625 


0.10983 


0.105719 


0.03125 


0.117483 


0.11255 


0.015625 


0.125715 


0.13276 


0.0078125 


0.134227 


0.140261 


0.00390625 


0.14287365 


0.147476 



TABLE III: Self-energy for the auxiliary field, evaluated at p = for L = oo and mL 



propagators will be involved in subdiagrams. Cancellations that in the continuum make use 
of equalities of counterterms will no longer hold. Thus we expect that nonsupersymmetric 
renormalization of the on-shell formulation first appears at two loops. 



E. Locality 

As stated above, we have computed Z^ — 1 = limp^o ^^{p)/p^- The fact that this has a 
finite limit demonstrates that at small p, the self-energy is analytic in p"^. Thus we find that 



and there is no evidence of nonlocality in the scalar self-energy. 



(4.18) 



V. SUPERCURRENT, MIXING AND RENORMALIZATION 



For a general superpotential W{(j)), the supercurrent is 



dW fdW 



(5.1) 



and in our case dW/dip = m(f) + g(f)'^. Because of the supersymmetry breaking on the lattice, 
this will mix with other operators in the same symmetry channel. For example, one has at 
the same naive dimension the operator 



T'^ = 9^0*P_x + d''(j)P+x 



(5.2) 



If the lattice action (Eq. (I2.22p in the massive case and Eq. (I2.23P in the massless case) is 
fine-tuned, then in the long distance effective theory there will be a supercurrent that is 

19 



conserved in the continuum limit. Thus one way to detect supersymmetry is to consider 
hnear combinations of bare lattice operators and extract the one that has vanishing four- 
divergence in the supersymmetric limit, modulo contact terms. We will study that approach 
in Section IV Dl below. Before doing so, we briefly consider a naive discretization of the 
continuum supercurrent (15.11) . 



A. Naive lattice supercurrent 

We formulate this lattice supercurrent in terms of tilded fields: 



dW fdWY 

D.^j^p+x + D,^*rp.x + -^I'P-x + Kr ^'^+* 

ocp \0(j) J 

—^ = m^ + g4>^ (5.3) 

acj) 

For instance the correlation functions 

(5,(x)0(y)x(O)C), (^,(a;)0*(y)x(O)C) (5.4) 

give rise to tree-level diagrams from the quadratic terms in S^. It is straightforward to 
obtain at this order 

(S^(x)0(l/)X(O)C) = V2^,^,P^S{x)D^^jG{y - x) + V2m*^,P^S{x)G{y - x) (5.5) 

where S{x) = {x{x)x{0)C) is the free theory fermion propagator and G{x) = (0(x)0*(O)) 
is the free theory scalar propagator. This will have the correct continuum limit 9^ 5^ = 0, 
modulo contact terms, since no divergences appear at tree-level. This can be explicitly 
checked by differentiating the expression above and using the equations satisfied by the free 
propagators: 

DiS{x) = -{mP+ + m*P_)S{x) + 0{a), x^O 

DiDiG{x) = \m\'^G{x) + 0{a), x^O (5.6) 

The variation of the action under the lattice version of the supersymmetry transformation 
(12A4D and (12J[5|1 is 

6S = V2a^J2x^'^ \9P+{'^<i>Di^ - Di^"^) + g*P-{2(t>*Di4)* - Di4)*'^)\ e (5.7) 

Note that this is 0{a), since 

lim Y^ X^CP+{2^Di4) - Di0^) = (5.8) 



X 



Also note the presence of g in (15. 7p . In order to see violation of the supersymmetric identity 
dfiS^ = in the continuum limit, diagrams involving the coupling g must be included. Loop 
diagrams are required, in order to get the diverenges that cancel the 0{a) factor coming 
from (15. 7p in the continuum limit. 
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B. Exactly conserved supercurrent? 

Because of (15. 7p . one can ask whether in the g = hniit there is an exactly conserved su- 
percurrent; by a proper latticization of the continuum supercurrent and the four-divergence, 
we would hope to find a discretized version of d^S^ = that holds at finite lattice spacing. 
In order to search for this hypothetical supercurrent, we follow the usual method and sup- 
pose that the parameter of the transformation e depends on spacetime position, e = e(x). 
Variation of the action with site-dependent e can generally be written in the form: 

A5 = -a^^e^(a;)C6(a;) (5.9) 

X 

Thus in the form of the action with all fields, Eq. (12.101) . 

1 SS 
& = 



a^ (5(e^C) 

= -V2{D,(PDiP^x + A'0^-X + DkPD^P+x + D2(l)DiP+x 

+D,(l)*D,P+x + Dl(t>*P+X + D^<P*D2P-X + D2(j)*D,P_x 

+D2FP+X - FD2P+X + D2F*P_x - F*D2P-x) 

-V2m{D^4)P+X + 4>DiP+x) - V2m*{D^4)*P_x + 4>*DiP-X) 

-V2g{24>D,4>P+x + 4>'DiP+x) - V2g* {24^* D.^ P^X + ^'DiP-x) 

+ -V2(Di<l>P+X + <^DiP+X + Di<t>*P_X + <^*DiP^X) (5.10) 

a 

On the other hand, when we just work with tilded fields, Eq. (I2.17p . 

6 = -v^(Di0^P_x - FpP_x + D,4>*pP+x - F*pP+x) 

+V2{-D,4>P.x + -D24>*P+x) 
a a 

-V2{D,P_x{l - ^D^r'F + DiP+x(l - ^D,)-'F*) 

-v^m(Di0P+X + 4>P>iP+x) - V2m*{D^4)*P^x + 4>*DiP-X) 
-V2g{24>D,4>P+X + 4>'DiP+x) - V2g*{24>*Di4>*P,x + 4>*'D,P^x) (5.11) 

Finally, for the quadratic terms there exist lattice identities analogous to integration by 
parts, such as: 

J2 Di<I^DiP-X = -J2 ^DlP-X. Y. D2ct>D,P+x = Y. ^DiD-^P+x (5-12) 

X XXX 

These sorts of identities give us, for either form of 6, 

^a^e = -V2J2(^'{9{'^^D^4>P+x + ^'DlP+x) 

X X ^ 

g* (20*Di0*P_x + 0*'i?iP-x) } (5.13) 
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In the continuum limit this quantity also vanishes; however, on the lattice we cannot inte- 
grate this cubic quantity by parts and the expression inside braces gives an 0{a) deviation 
from zero. 

In conclusion we do not find an exactly conserved supercurrent at g = but we do find 
an expression for which J2x © = 0, behaving as if © ~ •9^5'^, for g = 0. The quantity & 
will play an important role in the broken supersymmetry Ward-Takahashi identities of the 
lattice theory, which we consider next. 

C. Broken Ward-Takahashi identities 

As we have just seen, the variation of the action under an x-dependent spinor parameter 
ea{x) cannot be written as a simple product of e and a finite difference operator acting on an 
expression like (15. 3p . Otherwise f l5.13p would vanish. Rather, it takes a more general form, 
which we have denoted above as (15. 9p . Thus & will appear in the broken Ward-Takahashi 
identities that we are about the derive. Parenthetically, in the classical continuum limit, 
with the discretization (15. 3 p of the supercurrent, 

lim(6„(x) - d^S^,o.{x)) = (5.14) 

because we know that each term should go over to the continuum expression d^S^^a{,x). 
In addition to (15. 9p . we also must consider the variation of the source terms, 

(5.15) 



S^^^ = -a^ J2 \fCP+x + f^CP_x + ■J*4> + J4>* + k*F + kF* \ 



For this purpose we define 5^0 = {e^C)aAa4>, etc. The supersymmetry transformations of 
the tilded fields are given in (I2.16p . Thus 

A,X/3 = -V2{P4Di4> + F)C)p^ - V2{P4Di4>* + F*)C)p^, 

A J = V2(P+x)a, A J* = V2{P^x)a 

A„F = V2{D^P+x)a. A„F* = V2{D,P_x)a (5.16) 



Then we have the identity 



&a{x) 



A^Z{J;x) = = f[d4>d4>*dxdFdF*]e-^-^ 

+ {fCP+A^x){x) + {f^CP^A^x){x) 

+{rAj){x) + {JAj*){x) + {k*A^F){x) + {kA^F*){x) 



(5.17) 



Here Jj collectively denotes all the sources and S is the action (I2.17p . Since the identity 
holds for all values of J', derivatives of AZ{J'; x) with respect to the sources also vanish. 
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Thus it is the generating function of lattice correlation function identities associated with 
the supersymmetry transformations (12.161) . We refer to these as the broken supersymmetry 
Ward-Takahashi identities. 

We could measure a lattice transcription of the Ward-Takahashi identities that vanish 
in the continuum theory. Their deviation from zero gives a measure of the amount of 
supersymmetry violation in the latticized theory. A simple example is the following. 

= 



+ T7^¥7^^T-^]^c.Z{J,x 



= {-a^&^{x)xp{y) - V2{{D,4> + F)C)p^{x)5{x, y)) (5.18) 

This leads us to the identity 

V2Y,ci\F{x))C^p = -Y,a\<B^{x)xp{y)) (5.19) 

X x,y 

In the continuum theory the right-hand side vanishes because (3 — t- d^S^. On the lattice we 
have instead the simplified expression (I5.13p . which is an 0{a) lattice artifact. Divergent 
behavior in {&a{x)xii{y)) can lead to a result which does not vanish in the continuum limit. 
We can measure the violation of the continuum supersymmetry Ward-Takahashi identity 
j d'^x{F{x)) = by computing either side of (I5.19P in our numerical simulations. This is 
done for the left-hand side in Section IVIIDI below. 

Consider again the most general lattice action consistent with symmetries, Eq. (12.22^ . 
By tuning the parameters in this bare action we expect to obtain the supersymmetric limit 
for the long distance effective theory. On the other hand, once this fine-tuning is performed 
the bare lattice action (I2.22p will not be invariant with respect to the supersymmetry trans- 
formation of the bare fields: 

5S = -a\^C^&^Q (5.20) 

Thus we continue to have modified identities, (15.17p . which do not look qualitatively different 
from those away from the supersymmetric point in parameter space. What is different is 
that in the long distance effective theory there is a conserved supercurrent, d^S^ = as an 
operator relation in the continuum limit. To probe for supersymmetry, we must construct 
this supercurrent, built out of the appropriate set of bare field operators, Eq. (I5.2ip below. 

D. Supercurrent formulation 

We write down all operators that have the index structure of S'^a(x); leading repre- 
sentatives are given in Table HVl We denote these as 0]^l^ where n/2 is the engineering 
dimension. A linear combination of these is the long distance effective supercurrent at lattice 
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spacing a: 

s.M= E E^r^(«)«^"^'^^'^i':.(?(^) (5-21) 

n=3,5,7,... j 

On dimensional grounds, by" (a) must be a dimensionless quantity, and is therefore a 
function of the renormahzed parameters of the long distance theory, rrira and Qr, in the 
infinite volume limit. At finite volume it may also depend on a/L. We know that 

\imd^S^a.{x)=0 (5.22) 

a— >0 

as an operator relation, provided the relevant and marginal counterterms in the lattice action 
fl2.22p are tuned properly. If the fermion mass is set to zero, then we only need to tune the 
lattice action fl2.23p . In practice, we will add a fermion mass term to f l2.23p and extrapolate 
to vanishing bare fermion mass. 

In numerical computations it is convenient to take the spatial transform so that we instead 
work with 



a^ 



Y, df^S^o^it, x) = a^Yl dtSoait, x) = dtQo.it). (5.23) 



Since S^q, is fermionic, an odd number of x fields must appear in nonvanishing correlation 
functions with dtQa{t)- For this reason, correlation functions should involve operators with 
dimension of at least 3/2. In what follows we will exclusively consider two-point functions 
with fermionic operators of dimension n/2 which we denote O^ , where n = 3, 5, 7 etc. In 
order to identify the supercurrent nonperturbatively, what one really wants to study is the 
family of correlation functions 

M^'")(t) = {dtOtif{t)Ot^'\0, 0)) (5.24) 

where 

<fW = «'E<f(^'^) (5-25) 

Truncating 71,171 < A^max? "we want to tune the parameters of the bare action (I2.23P such 
that MlJ^ it) has a null space in the collective column index B = {j,n) for each i,m and 
t, which together form a collective row index A = {i,m,t): 

lim MABbB = (5.26) 

for all A = 1, ... , Nc, with b nontrivial. 

Of course at finite lattice spacing, no amount of fine-tuning will restore supersymmetry. 
Thus what we seek is not (I5.26p . but instead MAsbs = 0{a). However we need to find bs 
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dimension 


operators 


3/2 


7m^±X 


5/2 


lf.(pP±X, 7f,<P*P±X, P±^^,X, cr^uP±dyX 


7/2 


l^<t^''P±X. 7f.r^P±X, 7M''P±X, d^cbP±x, d^rP±X, <t^P±d^X. ^^i^^X, ^^±9^X, 
<y iiudu(t>P±X, (^t^udu(l)*P±X, (po-fj.uduP±X, (t>*<^iiiyduP±X 



TABLE IV: Operators that can potentially mix with the supercurrent, up to 0{a) supressed higher 
dimensional operators. 



with the constraint ||6|| = 1 in order to make the tuning independent of the normalization 
of this vector. We thus seek to minimize 



^ = Y,\J2'''-^'''">\'- I^I^bI'^I (5-27) 



A B B 

with respect to b. This is repeated at various 1712 and Ai until the minimum is found with 
respect to these parameters for fixed mi. Finally, we make an extrapolation to mi = 0, and 
identify the fine-tuned pair m,2 and Ai where F approaches its minimum value. 

We have taken ||6|| = 1 in our considerations so far. But we know that is not the whole 
story in the a — > limit. The operator S^ will undoubtedly need to be renormalized, and that 
will involve a divergent factor: S*?^^"' = ZsS^. This can be determined through a position 
space scheme. We demand that the free theory results are matched at some distance scale 
r: 

Zl{S,{x)SM)\.\=r = (^.(^)^^(0))free,|.|=. (5-28) 

Similarly, for the operators that appear in the correlation functions with S^ we must impose 

{Zt"'f{0t''\x)0\-''\(^))y^^^, = {Ot"\x)0\-''\(^)),^^^^^^^^, (5.29) 

Then the quantity that should be minimized, and which must vanish in the continuum limit 
is: 

F=Y, I ZsZt"^ J2 t^"^ (^*^K^ (t)Of/') (0, 0)) I ' (5.30) 

t,m,i j,n 



E. Continuum limit 

Since the model is a |0|^ theory coupled to fermions through a Yukawa interaction, it is 
expected to be trivial. The only continuum limit is therefore the free theory. We can work 
at arbitrarily small but finite lattice spacing a. Having tuned to the supersymmetric point in 
parameter space there is only one renormalized mass parameter rrir in the theory, and so in 
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lattice units we will measure rrira from the exponential decay of correlation functions of the 
elementary fields. Due to the U{1)r symmetry, this mass will be proportional to the bare 
fermion mass, which in lattice units is mia. It follows that we take a smaller by reducing 
the size of the lattice parameter mia. The triviality of the theory is then the statement that 
if mia = 0, the long distance effective coupling gr = 0, independently of the bare coupling 
Hi (recall that Ai is fine-tuned to obtain the supersymmetric limit). 

VI. FINE-TUNING WITH MULTICANONICAL REWEIGHTING 



Multicanonical methods 20|, |30|, l3l|] combined with "Ferrenberg-Swendsen reweighting" 
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- |23| have proven to be a powerful tool for maximizing the usefulness of Monte Carlo 
simulations over a range of parameter space. We refer to this combination of techniques as 
multicanonical reweighting (MCRW). For instance, MCRW was applied in a study comparing 



SU{2) and 5*0(3) = SU{2)/Z2 lattice gauge theories [32|, |33|. It was found to dramatically 
flatten the distributions with respect to three parameters, twists on gauge fields at the 
spatial boundaries. Another successful application of MCRW consists of lattice results for 



the electroweak phase transition 34. |35|. 



We will begin by describing MCRW generally for a theory of a real scalar field 0, followed 
by a presentation of how it would be applied to the lattice four-dimensional Wess-Zumino 
model that we are studying. 

A. Preliminaries 

Suppose we perform a Monte Carlo simulation at one value mo of the scalar mass m, so 
that the configurations sample the distribution determined by the action 

S{m = mo) = S{m = 0) + - / d'^x ■ml(f)^{x). (6.1) 



Following the "Ferrenberg-Swendsen reweighting" method 2lN23| one can use the follow- 
ing reweighting identity to compute the expectation value {O)^ of an operator O for the 
distribution with a mass m: 

^ ^"^ (e-^^^-^^Lo ^ Ec.n„)exp[-|(m2-mg)/#a:0y " ^"^ 

In the first equality (■ ■ ■),^^ is the expectation value with respect to the canonical distribution 
corresponding to (16.11) and 

AS{m) = -{m^- ml) f d^x 0^ (6.3) 

is the shift in the action when the mass is changed. In the second, j d^x 0^ and Oc are 
the mass term operator evaluated on configuration C and J2ceFin) ^^ ^^^ ^^^ °^^^ ^^^ 
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distribution F{n) of n configurations generated in the Monte Carlo simulation. These of 
course provide a finite ensemble that approximates the canonical distribution corresponding 
to ( I6.ip . The advantage of this approach is that one need only perform a single simula- 
tion at mass mo, storing the values of J d^x (pQ and Oc for each C, and then (C)^ can 
be computed for a swath of the parameter space m without having to perform any new 
simulations. Typically the time for this "offline" calculation is negligible compared to that 
of the simulation. 

Unfortunately, the regime of utility for this technique is limited by the overlap problem, 
in a way that often worsens exponentially in the spacetime volume. For instance, suppose 
the theory (16. ip has a quartic interaction and a critical mass-squared rn^ such that for 
m^ < ml there is spontaneous symmetry breaking. If we simulate with m^ > ml then the 
field is exponentially weighted toward j d^x 0^ ~ 0. Now suppose we attempt to reweight 
to rn^ < ml- In that case —{rn^ — rn^) > so that the exponential weight factor in (16. 2p 
is minimal at J rf^x 0^ ?^ 0. The ensemble that is generated in the Monte Carlo simulation 
will have exponentially few configuration in the regime where j d'^x cfP' is far from zero 
and e"'^'^'^™') is large. Because we will have very few representatives of configurations with 
the largest weight e"^'^*^™-', and most members of the ensemble have very small weight, 
fluctuations will be large and huge samples are required in order to have acceptable errors. 
The mismatch of the distributions gets worse as the number of lattice sites increases, because 
the exponent is extensive (i.e., scales like the spacetime volume L^ x T). 



As a concrete example. Fig. 4 of 32| shows that in the range of a three-dimensional 
parameter space the ordinary canonical Monte Carlo density of states varies by 14 orders of 
magnitude. This is for an 8^ x 4 lattice, which is still relatively small. The problem will get 
exponentially worse on larger lattices. 



In a number of contexts the technique of multicanonical reweighting 20|, l30|, ImI] has been 
found to ameliorate the overlap problem. One replaces S with 

Smcrw = S + W[0u02,...], (6.4) 

where W[Oi, O2, ■ ■ ■] is a carefully engineered function of some small set of observables. For 
instance in the model that we are studying, W will be a function of the mass term and the 
quartic term. 

0, = a'J2 \<f>\"^ ^2 = a^E |0r- (6.5) 

X X 

The (reweighted) expectation value of an observable in the distribution corresponding to 
Smcrw is: 

<^) = ^^ ,w[oo,oC] ■ (6-6) 

Since the e^ factor in (16. 6 p just cancels the e~^ Boltzmann factor coming from (16. 4p . 
one might wonder why it is introduced in the first place. The point is that the additional 
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Boltzmann factor e~^ in effect produces a weighted average over a continuum of canonical 
ensembles (hence the appelation "multicanonical" ) such that there is a good overlap with the 
distribution that one is reweighting to. The challenge is to design a W such that sampling 
is flattened over the range of observables one is interested in. 

We return to Fig. 4 of [32l]. It shows that the multicanonical Monte Carlo sampling 
distribution is fiat in the range of three-dimensional parameter space between the peaks, 
where the ordinary canonical Monte Carlo distribution varies by twelve orders of magnitude. 
The reweighting function W was represented by a numerical table, composed of the inverse 
density of states with respect to the tuned parameters. This is for an 8^ x 4 lattice, which 
is still relatively small, and it indicates that (9(10^^) more samples would be required in the 
canonical Monte Carlo approach in order to scan a comparable range of parameter space by 
ordinary Ferrenberg-Swendsen reweighting techniques. 



As another example, in studying first order phase transitions (e.g., 3J]), one chooses 
Oi to be the order parameter of the transition; in a model with a scalar field, typically 
Oi = J d^x 0^. One tunes iy[Oi] to cancel the nonperturbative effective potential for this 
operator, so that the Monte Carlo simulation samples evenly in Oi. This enhances statistics 
for configurations intermediate between the phases. In the mass scan example of Eq. (16. 2p . 
one has 

EceFjn) Oc exp {W[J d'x 0^] ~ U^l - mj) J d'x 0^) 

In this way, wherever the exponential in (16. 7p happens to be at its maximum, a large number 
of configurations will be generated, due to the fiat distribution with respect to J d'^x 0^. 
Two approaches exist for engineering a good function W. 

(1) One can employ a bootstrap method that iterates between Monte Carlo simulation 
and adjustments to W. For instance a numerical tabulation of density of states p may 



be obtained from a canonical simulation, as was done in [32|, |33[. Schematically, one 
obtains a histogram estimate of p(Ci) for an operator value range Oi range Oi,min < 
Oi < Oi,max- This provides an initial version of W, through W = l/p{Oi). If 
necessary, the process can be repeated to refine the table. 

(2) Iterative or stochastic searches may be used to optimize W with respect to a prede- 
termined parameterization in a small volume. Performing this at two different small 
volumes then provides an extrapolation estimate for W in the next largest volume, 
which can then be refined through another search. 

B. Application 

If we work at rrii ^ 0, the action that we must tune nonperturbatively is (I2.22p . If 
extrapolate to the mi — > limit, t/(l)ij symmetry (12.131) can be imposed and we must 
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fine-tune fl2.23p . We fix yi, so that it determines the coupling strength that we study. All 
other parameters, m^ and Ai, are associated with bosonic terms in the action, and can be 
tuned offline in the way that was just described above. 

VII. SIMULATION 

A. PfafRan phase 

The integration over lattice fermions yields 

Pf (CM), M = p + mP+ + m*P^+g4)P++g*4)*P_ (7.1) 

The transformation 

X(a;°,x)^7°x(a;°,-a;), 0(a;°,a;) ^ 0*(x°, -a;), m^m*, g ^ g* (7.2) 

leaves x^CMx invariant. Its effect on the Pfaffian is then 

Pf(CM)(m,(?,0) = Pf (CM)(m*, (7*, 0*) = [Pf(CM)(m, (7, 0)]* (7.3) 

so that the Pfaffian is real. There is still the possibility of a sign problem. At weak coupling 
and \m\ not too small this is likely circumvented, since the spectrum of M is pushed away 
from zero and crossings are avoided. In that case we can work with the phase quenched 
measure, 

|Pf(CM)| =det(M^M)i/^ (7.4) 

B. Pseudofermions 

The Pfaffian is represented through an integration over bosonic fields, the pseudofermions 
?7, with action 

SpF = r]\M^M)-^'^T]. (7.5) 

We compute this using the rational approximation 

d 



(MtM)-^ - «or/ + $: -TTT^-jV {7.Q) 



i=l 



where the coefficients ctj, /3j are chosen to minimize errors over the range of eigenvalues of 
the operator M'^M. Since we will work at weak coupling in what follows, we are able to use 
the spectrum of the g = operator, MgMo, to determine this range. In (17. 6p . the integer d 
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Lattice 


CPU 


GPU (CUDA) 


GPU (Ours) 


CPU 


GPU (CUDA) 


GPU (Ours) 




single 


single 


single 


double 


double 


double 


8^ x32 


1.1 


6.9 


24 


0.94 


4.4 


11 


16^ X 32 


0.88 


14 


71 


0.69 


10 


35 


323 X 64 


0.11 


20 


— 


0.085 


10 


— 



TABLE V: Comparison of timing, Gflop/s, for fast fourier transform of the spinor field. Both 
single and double precision results are given. 



is the degree of the approximation. Evaluation of (17. 6p requires us to solve the linear algebra 
problem 



{M^M + /3,)X, = r] 



(7.7) 



for each degree i = 1, . . . ,d. This is done for each i by conjugate gradient, which is an 
iterative solver designed for sparse linear systems. In computing the M'^M matrix multipli- 
cation, we encounter Mfin = font and a similar expression involving MK There is a trick for 
the computation of the "dslash" p; we fast Fourier tranform the in vector /j„ to momen- 
tum space, apply dslash here where it is diagonal (giving font in momentum space), then 
fast Fourier transform back to position space. This avoids an additional layer of rational 
approximation for computing products involving (A^AY^'^ in position space. 

In Table |V] we give a comparison of timing for fast fourier transform (FFT) of a lattice 
spinor field such as /j„ using different tools. Concerning hardware, the CPU runs are from 
an Intel Xeon Woodcrest 5160, whereas the GPU runs are from a Nvidia GTX 285. The 
CPU runs were performed using the four-dimensional Numerical Recipes code 36|]. The first 
version of GPU code performs four one-dimensional transforms using the batched CUDA 
FFT available from Nvidia to get the four-dimensional transform. The second version of 
GPU code uses our own FFT for the spatial dimensions, which is currently only operational 
for L < 32, due to register constraints. As can be seen, it is quite a bit faster than CUDA 
FFT. In the temporal dimension T > 32 so the CUDA FFT is always used in that direction. 
A necessary reorganizing of the arrays after each of the four one-dimensional transforms 
is included in the timing. We do not know the origin of the drop in performance of the 
Numerical Recipes code on large lattices, though it may be related to non-cached memory 
accesses. We inserted a counter of floating point operations (flop) into the Numerical Recipes 
code and found that the number of operations was approximately given by 5 A^ log2 N. We 
used the estimate 5A^log2A^ for the number of operations for GPU code as well. What 
one sees is that even on small lattices the FFT runs an order of magnitude faster using the 
GPU. Since this is the crucial operation in applying the lattice Dirac operator, we anticipate 
simulation speeds of approximately 10 Gflop/s for the 16^ x 32 and larger lattices, working 
in double precision. In what follows we use our own FFT for the L < 32 lattices, but the 
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CUBA FFT for the L = 32 lattice. 

C. Preconditioning 

At weak couplings, we expect that preconditioning by the inverse of the free theory 
fermion matrix Mq will improve convergence of the conjugate gradient solver. We must 
formulate the problem in terms of a hermitian, positive definite matrix. For this purpose, 
we re-express the problem 

M^Mx = b (7.8) 

as follows: 

{M-^^M^){MMq^){Mox) = (Mo^^b) (7.9) 

This leads to the definitions 

M = MMq\ x = Mox, b = M~^^b (7.10) 

The problem in terms of these variables is M'^Mx = b. Once we obtain x from the conjugate 
gradient solver, we obtain the desired solution from x = M^'^x, where Mq^ can be computed 
analytically. It is just 

Mo-'= .^^°, „ , p^j,D, = il-^D,rD^ (7.11) 

We have solved the problem M'^Mx = b with random Gaussian b on lattices of various 
sizes, using our GPU code. Convergence is obtained with and without preconditioning. A 
comparison is given in Table IVII The speed-up factor from the preconditioning is quite 
large, and the inversion times are rapid enough that simulations on relatively large lattices 
are realistic. 

The pseudofermion rj is updated by the heatbath method, 

r]={M^MY^'^R (7.12) 

where i? is a random complex Gaussian spinor. The 1/8 power is obtained by rational 
approximation so that rather than the problem (17.81) . we must solve 

{M^M + /3i)x = b (7.13) 

for each degree i = 1, . . . ,d. In order to have high accuracy we have taken d = 20. Applying 
the preconditioning matrix Mq~^ then yields 

M^M + PiM^^^M^^^ x = b (7.14) 
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Lattice 


Precision 


NPC sees. 


NPC iters. 


Gflop/s 


PC sees. 


PC iters. 


Gflop/s 


speed-up 


8^ x32 


single 


1.3 


830 


14 


0.038 


8 


17 


34 


8^ X 32 


double 


4.9 


1600 


7.2 


0.13 


15 


8.1 


38 


16^ X 32 


single 


7.1 


870 


25 


0.20 


8 


31 


36 


16^ X 32 


double 


31 


1800 


12 


0.74 


15 


13 


42 


32^ X 64 


single 


420 


1600 


14 


6.8 


8 


15 


62 


32^ X 64 


double 


2200 


3900 


6.5 


25 


15 


7.0 


86 



TABLE VI: Timing benchmarks at m = 1, g = 1/5. PC indicates preconditioning whereas NPC is 
the inversion without preconditioning. The time in seconds and the number of iterations (iters.) 
is for convergence. The criterion for convergence is that the relative residual is less than 1 x 10^^ 
for single precision and less than 1 x 10^^^ for double precision. The speed-up is the ratio of NPC 
time to PC time. 



Lattice 


Sees. 


Gflop/s 


Iters. 


Sees. 


Gflop/s 


Iters. 




single 


single 


single 


double 


double 


double 


8^ x32 


0.93 


15 


360 


3.5 


7.5 


700 


163 X 32 


4.8 


28 


370 


20 


13 


730 


32^ X 64 


220 


14 


540 


900 


6.7 


1100 



TABLE VII: Pseudofermion heatbath timing benchmarks atm = l,g = 0.1, with degree 20 rational 
approximation of (|7.12p . Preconditioning is used up to /3j = 100. The number of iterations is the 
total over the 20 different degrees. 



We find that there is degradation of the speedup due to preconditioning for larger values of 
/3j. However for very large (3i the problem without preconditioning converges rapidly. Thus 
we choose a cutoff value, which turns out to be (3i = 100, above which preconditioning is 
not used. By this approach we are able to keep the number of iterations for each of the 
problems f l7.13p under 100 for single precision and under 200 for double precision. As Table 
IVIII shows, the average number of iterations is far less. For example, at single precision on 
the 8^ X 32 lattice, a total of 370 conjugate gradient iterations are needed over 20 degrees, 
for an average of 18.5 iterations for each problem fl7.13p . 

In the molecular dynamics part of the rational hybrid Monte Carlo (RHMC) algorithm, 
the fermion force terms also require the solution of equations of the form f l7.13p . These have 
a performance similar to what is shown in Table IVIII We find that approximately 50 steps 
with double precision are required in order to get good acceptance rates for the Metropolis 
step that is applied at the end of a molecular dynamics trajectory of length 0.5 simulation 
time units. We have found that single precision is inadequate for the molecular dynamics 
evolution. It could however be used for a mixed precision conjugate gradient, as has been 
promoted in 37 . 
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D. Ward-Takahashi identity 

In the continuum the simplest Ward-Takahashi identity is that 

(F) = (7.15) 

On the lattice this is modified due to the noninvariance of the lattice action. It is therefore 
of interest to measure (F) in our simulations, which as alluded to above use the RHMC 
algorithm. Working at ma = 0.5, gf = 0.1 we find that sampling over 3283 configurations, 
after 200 thermalization sweeps, 

^ E (^(^)) = [(2 ± 3) + ^{2 ± 3)] X 10^^ (7.16) 

The error estimates incorporate an autocorrelation time of approximately 15 sweeps, for this 
observable. The small value of {F) is consistent with the fact that the coupling is weak and 
at one- loop order in perturbation theory, nonrenormalization of the superpotential is intact. 
Thus we expect corrections of the scalar potential to start at Od^fl^) = (9(10~^), roughly 
consistent with the size of (F). Nonsupersymmetric corrections to the scalar potential are 
required in order to have (F) nonzero. 

VIII. CONCLUSIONS 

We have studied the locality of the Dirac operator that appears in the action 02. 171) with 
fermion matrix (I2.19p . It was seen in Fig. [1] that there is a long tail on this operator, so 
that the localization is less pronounced than that of the overlap operator. We found that 
the scalar self-energy at one-loop takes the form (I4.18p . so that there is no indication of 
nonlocality in the external momentum. Thus it may be that the long tail of the operator is 
harmless. 

We established an equivalence between the "untilded" formulation of this lattice theory 
and the "tilded" theory, at the level of correlation functions in the continuum limit, provided 
all divergences are logarithmic. 

We studied one-loop counterterms in the lattice theory. We have found agreement with 
previous authors that in the quantum continuum limit of this theory at one-loop, only 
wavefunction renormalization occurs. We have given numerical values and find that in 
agreement with what appears in the literature, Z^ = Z^ but that Zp is different, in the 
continuum limit at one-loop. 

We described the renormalization of the bare lattice supercurrent. An approach to de- 
termining the linear combination that is the symmetry current of the long distance effective 
theory was outlined. This is crucial for fine-tuning the bare lattice action in the reweighting 
approach that we are pursuing. We showed that the invariance of the free lattice action does 
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not have a simple interpretation in terms of a conserved supercurrent, due to the fact that 
the derivative operators in this formulation are not ultralocal. 

We have developed code to run on graphics processing units that are CUDA enabled, and 
have studied a strategy for handling the inverse square roots of matrices that occur in the 
theory through fast Fourier transform libraries within CUDA, as well as our own improved 
version which runs on smaller lattices. We have shown that preconditioning with the inverse 
of the free theory matrix yields a speed-up that is quite significant in the problems f l7.13p 
that are encountered in an RHMC sweep. 

In future work we will fine tune the bare lattice action f l2.23p through measurements of 
indentities involving d^S^, taking an approach that involves reweighting and multicanonical 
Monte Carlo techniques. We have already measured (F) and find it to be quite small at 
weak coupling. Thus the bare theory without any fine-tuning is a good starting point for 
the search over parameter space. An incremental search strategy, starting at weak coupling, 
moving gradually to stronger coupling, looks realistic. 
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